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The effects of intrinsic noise on stochastic delay systems is studied within an expansion in the 
inverse system size. We show that the stochastic nature of the underlying dynamics may induce 
oscillatory behaviour in parameter ranges where the deterministic system does not sustain cycles, 
and compute the power spectra of these stochastic oscillations analytically, in good agreement with 
simulations. The theory is developed in the context of a simple one-dimensional toy model, but 
is applicable more generally. Gene regulatory systems in particular often contain only a small 
number of molecules, leading to significant fluctuations in mRNA and protein concentrations. As 
an application we therefore study a minimalistic model of the expression levels of hesl mRNA and 
Hesl protein, representing the simple motif of an auto-inhibitory feedback loop and motivated by 
its relevance to somite segmentation. 
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I. INTRODUCTION 



Most processes in biology are intrinsically stochastic, due to the random fashion in which molecules interact. 
In order for a biochemical reaction to occur, for example, all reagents must be sufficiently close in space, and 
due to thermal or other types of stochasticity the execution of reaction is fundamentally a stochastic process. 
This type of randomness has, until recently, mostly been neglected in attempts to model biochemical reaction 
systems, and deterministic rate equations have often been used to describe the dynamics of such systems. Noise 
has here often been assumed to have only a minor effect on the dynamics, so that it could safely be ignored. 
The use of deterministic approaches implies the assumption of large, formally infinite system sizes, only in this 
limit can the law of large numbers be applied to show that the resulting mean-field equations give an accurate 
description of the dynamics of the system. Additionally, it is frequently assumed implicitly that the reactor 
^f-^ ■ in which the chemical dynamics takes place is well-mixed, so that spatial variation of concentrations of the 
interacting chemicals can be ignored. 

The reason for the popularity of such approaches undoubtedly rests in their relative mathematical simplicity: 
while the methods with which to analyse sets of non-linear deterministic differential equation are fully developed 
(see e.g. [l| or similar textbooks), a theory for the corresponding stochastic systems is far less advanced. If 
the number of reacting molecules in the system is small, the stochastic effects can no longer be ignored. An 
important example are mRNA molecules in gene expression [1, [!, 0, Hj] , where only a small number of molecules 
is involved in the reaction dynamics. Unsurprisingly deviations from the mean-field dynamics are here to be 
expected, and stochastic rather than deterministic modeling approaches to such systems in molecular biology 
are appropriate. Only in recent years have analytical and more systematic studies of such systems been 
undertaken, and substantial differences between the behaviour of stochastic systems and their deterministic 
counterparts have been found in different model systems. In particular so-called demographic noise |6[ due to 
the discreteness of the dynamics may change the structure of the attractor of a given system fundamentally. 
References relevant for the present work can be found e.g. in [1, 0, H, 0, H, @ • 

Stochastic approaches to biochemical reaction systems typically start from a master equation describing the 
microscopic dynamics, the deterministic mean-field dynamics can then formally be derived to lowest order 
within a van Kampen expansion in the inverse system size [13]. Taking into account next-to- leading order 
finite-size corrections can alter the dynamics considerably. Predator-prey systems with a fixed point on the 
deterministic level can for example be seen to exhibit coherent sustained oscillations at finite sizes Q. The 
spectrum of these cycles can be obtained to striking precision within the system-size expansion. Similar 
oscillations have been found in a variety of other systems, including models of epidemiology, opinion dynamics 



* Electronic address: tobias.gallaOmanchester. ac.uk 



2 



and biochemical reaction networks [8|, |9| . 

In addition to the discreteness of the dynamics and the resulting intrinsic stochasticity, processes in gene 
regulatory systems are typically subject to considerable delays induced by the underlying biochemical reactions. 
I.e. processes such as transcription and translation do not occur instantaneously, but generate their reaction 
products only well after the reaction has been triggered [1, H, @ . The aim of the present paper is therefore 
to extend the theoretical tools developed in [7|, la, |9( to the case of stochastic delay systems, and to use them 
to study a simple model of gene regulation. As we will see the dynamics of such systems may well exhibit 
stochastic coherent oscillations at finite system size in ranges of the reaction rates where the deterministic, 
infinite delay system has no cycleSjbut approaches a fixed point instead. Such oscillations in delay stochastic 
systems have been reported in 0, 0, IE] , but to our knowledge a theoretical computation of correlations and 
power spectra of these cycles within a systematic expansion in the inverse system size has not been attempted 
in the context of such models. Theoretical approaches based on generating functions have however been 
discussed in Q . In this paper we will follow a complementary approach, and in particular we will describe how 
the method of the system-size expansion applies to delay systems, and how a linear delay Langevin equation 
can be derived for fluctuations about the trajectory of the deterministic mean-field system. From this Langevin 
equation the power spectra of these stochastic cycles can then be obtained analytically. Our analysis therefore 
offers a theoretical characterization of results from simulations reported e.g. in [4|, and an alternative to the 
analytical approaches discussed in 



II. TOY MODEL 



A. Definition and deterministic description 

In order to develop the formalism we start with a simple model of delay stochastic processes, and consider 
a system in which there is only one reacting substance X. The dynamics are given by the following reactions 

A — > A + X, (1) 
B + X — > B, (2) 
X + C =>■ C, (3) 

note that neither reaction affects the number molecules of substances A, B and C in the system, so that these 
reactants are mere 'dummy' variables, and their only role is to set the relative rates with which the three 
reactions occur. While the first two reactions are assumed to occur instantaneously, the third reaction involves 
a delay, we indicate this by the double arrow in Eq. ([3]). I.e. if a reaction between a molecule of type C and 
a molecule of type X is triggered at time t, then one molecule of type X is removed from the system at a 
later time t + r (provided there is at least one X-molecule in the system at this later time) . For simplicity 
we will assume that r is a constant delay period, but variable delay times, drawn e.g. from some probability 
distribution can be in principle considered as well 0. The model in this setup has previously been discussed 
and studied within an alternative approach in 

On the mean-field level the concentration of X-molecules in the system is described by the delay differential 
equation 

x(t) = a - bx(t) — exit - t), (4) 

where a, b and c are non-negative rate constants related to the (constant) number of molecules of types A, B 
and C in the system respectively (further details will be discussed below). 

Eq. (j4]) is linear and its asymptotic behaviour can be computed straightforwardly. In particular a linear 
stability analysis has been carried out in p|, and a phase diagram was obtained in terms of the parameters 
a, b, c and r. The only fixed point of Eq. ^ is x* — a/(b + c), and at fixed a it is found to be unstable at 
large rc or small rb respectively. In such circumstances oscillations grow indefinitely. Below a the line marking 
the Hopf bifurcation in the (rb, rc) plane, the fixed-point is stable [5(. We will focus on this regime in the 
following. 



B. Stochastic dynamics, van Kampen expansion and spectrum of fluctuations 



In order to model the dynamics on the microscopic level, let us assume the reactor in which the various 
reactions take place contains afl molecules of type A, &f2 molecules of type B, cil molecules of type C, and 
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n particles of type X. Since the number of A, B and C molecules remains unchanged by the reactions given 
above, the only dynamical variable in the system is n(t) (which will be of order f2 as are the numbers of the 
other particles in the system). The first reaction, Eq. (fTJ), then occurs with a rate aO, the second reaction with 
rate nb and the third with rate nc, and the resulting stochastic process is described by the master equation [5] 



—P(n,i) = ail(l 
at 



- l)P(n, t) + b(E - l)[nP{n,t)] 



+c V m(E- l)[Q(n)P(n,t;m,t - t)] 



m=0 



(5) 



for the probability P(n, t) of finding the system in state n at time t. E is here an operator acting on a function 
of n via E/(n) = f(n + 1), and not to be confused with an expectation value of some kind. E _1 stands for 
the inverse operation. 0(n) is the step function, i.e. Q(n > 0) = 1 and Q(n = 0) = which ensures that the 
delayed removal of X-molecules only occurs provided there is at least one molecule present in the system at 
the time at which the removal is due to take place. Note that Eq. (O is not closed on the level of one-time 
quantities, as P(n, t;m,t — r) describes the joint probability distribution of finding n X-molecules at time t, 
and m X-molecules at time t — r. 

Following (lOj and anticipating that n will be of order Q with fluctuations of order Q 1 / 2 we now introduce 
continuous degrees of freedom, and write 



n(t) 



x(t) 



Q 1 / 2 ' 



(G) 



The above master equation for P(n, t) can then be written in terms of the distribution II(£, t), and within an 
expansion in powers of f2 -1 / 2 we have similar to 10] 



W,*)-^^M : 



<9£ 



2de 



2de 



(nx(t) + n^ 2 ^u^,t) 



2de 



(7) 



where we introduce 77 by writing n(t—T)/Q = x(t — t)+t]/{1, 1 / 2 , and where we have ignored higher-order terms. 
Anticipating that we will take the limit of large systems eventually and that trajectories at which n(t) = 
at any time will not contribute in this limit we have ignored the factor O(n) in the last term of ([5]). This 
is common in the context of the van Kampen expansion, which is usually unable to capture features such as 
absorbing states at the boundaries of configuration space. Collecting terms of order fi 1 / 2 in Eq. ([7]) one finds 



n^x = -an^+bxm^+cxit-r)^ 2 ^ 



(8) 



in the lowest order of the van Kampen expansion. Wh have here written II as a shorthand for II(£, t) and used 
the identity II(£, t) = J drj II(£, t; r\, t — r). From (j8]) one has 



x{t) — a — bx(t) — cx{t — r), 



(9) 



i.e. one recovers the above deterministic equation ((3]). To next-leading order (i.e. collecting O(fl ) terms) one 
finds 



1 d 2 



hx(t)^U(U) 



(10) 



Integrating out rj in the last term one then has 
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FIG. 1: (Colour on-line). Dynamics of the toy model at a = 100, b = 4.1, c = 4.,r = 2. The dark line with decaying 
oscillations in the main panel shows the behaviour of the concentration of X molecules in the deterministic system, 
Eq. ([9| [111 ]. The noisy line with persistent oscillations represents one simulation run of the stochastic dynamics at 
Q = 100. The inset shows a zoom at large times in the equilibrated regime, the horizontal line is the mean-field fixed 
point, the stochastic system shows persistent cycles. 



At asymptotic times t the mean-field trajectory approaches its fixed point (our analysis is restricted to the 
stable phase, for similar studies in non-delay systems with a limit cycle see [HI). We therefore replace x(t) 
and x(t — r) in Eq. (fTT|) by x* = a/(b + c). With this substitution Eq. (fTTj) then describes a delayed Langevin 
dynamics of the form 

i = -b£(t)-c£(t-T) + <;(t), (12) 
where C(i) is Gaussian white noise of zero mean and with variance 

(C(*)C(0) = ( a + bx * + cx *) S ( t - *')■ (13) 

See e.g. Frank et al. [Hj] for further details on Fokker-Planck descriptions of delay systems. Eq. Q12j) is linear 
and can be solved in Fourier space (similar approaches to delay Langevin equations have been discussed in 
[HI]). In particular one has 

[icj + b + ce- iUT ]l{u) = CM, (14) 

where £(u>) and C(^) indicate the Fourier transforms of £(i) and ((t) respectively. From Eq. (fT4| one directly 
reads off the power spectrum of £ (t) and finds 

p(u) = (ichi 2 ) 

a + bx* + cx* 2a 



[b + ccos(wt)] 2 + [lu — csin(wr)] 2 [b + ccos(wt)] 2 + [u — csin(o;T)] : 



(15) 

We have here used the relation (^((u;)((uj')j = (a + bx* + cx*)8(u + to'), where (• • •) denotes an average over 

the stochastic process described by the Fokker-Planck equation (jTTJ) , or equivalently over realizations of the 
Langevin equation (fT2|) . 



C. Test against simulations 

Microscopic simulations of the processes defined by the reactions (QJ3]) can be carried out using the algorithm 
originally proposed by Gillespie [l5|], suitably modified to take into account the delayed reactions. Details of 
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such modified Gillespie schemes have been discussed for example in [161 ]. but for completeness we re-iterate 
them here. Essentially the simulations follow that of the classic Gillespie algorithm and whenever a 
delayed reaction is triggered it is added to a list of delay reactions to be executed at a later time (r units of 
time after the reaction is initiated). This list is constantly updated, and delay reactions are executed (and 
removed from the list) in a manner consistent with the probabilistic description in terms of the above master 
equation. Specifically the simulations of our toy model dynamics proceed according to the following algorithm: 

1. Initialise. Set model parameters a,6, c and the system size fi. Set the initial number n of molecules X, 
and set t = 0. Set list of scheduled delay reaction to an empty list. 

2. Calculate the propensity functions a\ — of2, a,2 — bn, <j(3 = cn. 

3. Compute ao = a>\ + a 2 + &3- 

4. Generate an independent random number r from a uniform distribution over (0,1], and set A = 
-ln(r)/a - 

5. If there is a delayed reaction scheduled to occur during the interval [t, t + A) then 

(a) Identify next delayed reaction scheduled, and, provided n > execute it, i.e. reduce n by one. If 
n = before the reaction, then do not execute the update (otherwise n would go negative). In 
either case remove the reaction from the list of scheduled reactions. 

(b) Update t to the time for which this reaction was scheduled. 

(c) Go to 2. 

6. If there is no delayed reaction scheduled for the interval [t, t + A) then 

(a) Generate an independent random number r from a uniform distribution over (0,1], and find \x G 
{1,2,3} such that 



(b) If /i = 1 or /i = 2 and then execute the corresponding reaction (not a delay reaction), and increment 
t by A. Go to 2. 



(after a suitable equilibration time) , where x* is the asymptotic fixed-point value of the deterministic dynamics 
given by Eq. @, i.e. x* — a/(b + c), or equivalently the long-time average of n(t)/Sl. From these time series 
£(i) a numerical measurement of the power spectrum P(u>) is obtained by subsequent Fourier transform, and 
finally results are averaged over a sufficiently large number of independent runs. 

Results of stochastic simulations of this system are shown in Figs. Q]and[2] The first figure depicts a single 
run of the stochastic system and shows that coherent oscillations are sustained in parameter regimes in which 
the deterministic equations approach a fixed point. The mechanism by which these oscillations are generated 
is the following: the deterministic system approaches its fixed point in an oscillatory manner (i.e. the Jacobian 
at the fixed point has eigenvalues with non-zero imaginary parts) , and the stochasticity of the finite system 
results in persistent perturbations away from this fixed point, so that both features together result in an 
overall oscillatory effect. This has been seen in a variety of different systems |7|, |9(, but it is worth pointing 
out that in non-delay systems a minimum of two dimensions is necessary to allow for a complex eigenvalue of 
the Jacobian. In delay systems one degree of freedom is sufficient [5(, so that even the one-dimensional toy 
model discussed in this section is able to produce demographic oscillations about the mean-field fixed point. 
Fig. [2] demonstrates that the analytically obtained spectrum, Eq. (|15p . agrees very well with simulations, we 
attribute the remaining small discrepancies to finite-size or equilibration effects. A similar figure was obtained 
by different methods (based on generating functions) in Q. 
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FIG. 2: (Colour on-line) Power spectrum P(ui) = ( |£(w)| 2 \ of the fluctuations about the mean-field fixed point of the 

toy model. Parameters a, b and c are as in Fig. [T] The line shows the analytically obtained spectrum of (|15[) . symbols 
represent results from simulations at Q = 200 (averaged over 341 samples, measurements start at t = 100 to allow for 
some equilibration period). 



III. SIMPLE MODEL OF GENE REGULATION 



The second system we will be studying is a simple model of gene regulation. We here chose a system that 
represents one of the most common motifs in gene regulatory networks, namely a model of a single gene- 
protein synthesis with negative delayed feedback [5(. Specifically we address a model previously discussed in 
[4j, describing the coupled time behaviour of the expression levels of so-called hesl messenger RNA (mRNA) 
and Hesl protein. Following the notation in the biology literature we will italicize and use lower case when 
referring to mRNA, and will use non-italicized font with the first letter in upper case when referring to the 
protein 0, [T3]. Hesl here is a Notch-signalling molecule, where the so-called Delta-Notch signalling process is 
a mechanism for cell-cell communication and underlies cell differentiation for example in vertebrates 0, [H, [lj| . 
Cycles with a time-period of approximately two hours has been reported for the concentrations of hesl mRNA 
and Hesl protein in mice [J, [13] ■ These oscillatory processes, also referred to as the somite-segmentation clock 
[lj|, are linked to the formation of somites, i.e. the emergence of blocks of cells which determine the future 
positions of skeletal muscles or vertebrae [H, [lji[ . Spatial segmentation in the body here can be understood 
as a reflexion of temporal oscillations in gene expression [3|, I la . 1 19L 1201 ] . The underlying molecular mechanism 
producing the oscillations of mRNA and protein are hence of great interest, and several mathematical models 
have been proposed, among them 0, H, 0, HH and [23 ]. See also 0, (23| for stochastic effects in models of gene 
regulation. 

We will not discuss the details of the biochemical mechanisms in this paper, but will only present a brief 
abstraction of the reactions necessary to define the mathematical model we will study here. Further details 
on modelling genetic circuits can be found in 24] or in similar textbooks. In essence the model describes the 
concentrations and interactions of two types of substances, hesl mRNA and Hesl protein, as illustrated in 
Fig. [21 mRNA molecules are produced by transcription of DNA. This involves several biochemical processes 
(e.g. elongation, splicing) which we will neglect in our description. The rate at which mRNA molecules are 
produced depends on the concentration of protein through a negative feedback mechanism. Transcription is 
here associated with a delay time r, so that it is the protein concentration at time t — r which affects the 
production rate of mRNA at time t. This will be explained further below. mRNA is also subject to degradation 
(i.e. removal from the system) at a constant rate /i TO . In a process subsequent to transcription hesl mRNA 
molecules are translated into Hesl protein (the mRNA molecule is not used up in this process). Translation 
may involve another delay, which for simplicity can be absorbed into the transcriptional delay Protein 
molecules finally are subject to a degradation process at rate \i v . Crucially, a negative feedback is induced 
by a repressatory effect of Hesl protein on the transcription process. Hesl dimers may bind to the relevant 
promoter regions in the DNA, and reduce the transcription rate at which mRNA is generated (following [4] 
dimerization is not discussed as an explicit step in our work, but see Q for models taking this into account). 
The transcriptional repressor Hesl thus negatively affects its own expression [13]. Mathematically this is 
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FIG. 3: (Colour on-line) Schematic illustration of the Hesl regulatory system. See [2ll ] for a similar picture. The 
first process (transcription) includes the elongation, splicing, processing and export from the nucleus of primary gene 
transcript. The synthesis of Hes 1 protein occurs by translation of hesl mRNA. Any translational delay is here absorbed 
into the transcriptional delay time r. The Hes 1 protein finally represses the transcription through the binding to the 
promoter. Both, the mRNA and the protein are subject to degradation. 



modelled by a transcription rate which depends on the concentration of protein via a decreasing function, as 
we will now explain. We will focus on the model proposed in [H, [HI]. In its deterministic form it is given by 
the differential equations 

j t M(t) = a m f(P(t-T))-» m M(t), (16) 

j t P{t) = a p M(t) - fi p P(t). (17) 

M(t) here labels the concentration of hesl mRNA, and Pit) that of the Hes 1 protein. /i m and \i p are 
degradation rates for the mRNA and the protein respectively, and a m is the mRNA transcription rate in the 
absence of protein expression (f(P) is still to be defined, but we will have /(0) = 1). a p is the translation rate. 
f(P) finally is a monotonically decreasing Hill function representing the suppression of mRNA production 
through the binding of Hes 1 protein dimers into the promotion region. In the model it takes the form [J, [24[ 

rem = TTmm (18) 

with h the so-called Hill coefficient. Po is the concentration of protein at which f(P = Pq) = 1/2. Eqs. (|16I17[) 

are the deterministic abstraction of an underlying microscopic stochastic model defined by the following four 
reactions [H 

M ^ 0, (19) 

P 0, (20) 

M M + P, (21) 

^4 M. (22) 

These dynamics may be described by a stochastic process for the numbers n m of mRNA molecules and n p of 
protein molecules in the system. For later convenience we will write n = (n m , n p ). The first two reactions here 
describe the degradation of mRNA and protein respectively. The third reaction captures the translation of 
mRNA into protein. The last reaction finally corresponds to the production of hesl mRNA via transcription. 
Note that DNA is not part of the dynamical model (its concentration is constant), which is why M appears 
to be produced out of the void in the fourth reaction. In absence of protein (n p = 0) this reaction occurs at a 
rate a m , but is suppressed by the Hesl protein, and in total the rate at which mRNA is produced at time t is 
hence Qt m f(n p (t — r)/0), where fl is a measure of the system size. The time lag r in the argument of / models 
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FIG. 4: Time series of the concentrations of mRNA and protein concentrations respectively. Solid curves, decaying 
towards a fixed point, are from a numerical integration of the deterministic dynamics, Eqs. (|16I17|I . Curves with 
persistent oscillations represent a single run of the stochastic dynamics at Q = 1000. Model parameters are as some 
of the examples in Q: Po = 10,h = 4.1, r = 18.7, a m = a p — l,fi m = A*j> = 0.03. Units of a p ,a m and of /i m ,Mp are 
min" 1 , r is measured in minutes [3]. 



the delayed repression of hesl mRNA production by the protein. The rate of production of mRNA at time t 
is therefore negatively regulated by the concentration of protein at time t — r. A typical run of the stochastic 
system is shown and compared with the deterministic system in Fig. [J] Model parameters are chosen as in 
[J|. As seen in the figure, the deterministic system approaches a stable fixed point asymptotically, with a 
complex eigenvalue as indicated by the decaying modulations. The stochastic system at finite size remains 
in an oscillatory state as previously observed in 0]. We will now proceed to characterize these oscillations 
analytically, applying the formalism developed in the previous section. 
The master equation describing the processes (|19H22j) then takes the form 

^P(n;t) = fx m (E M - l)[n m P(n;t)}+ fi p (E P - l)[n p P(n; t)] + a p (E p l - l)[n m P(n; t)] 

+a m nJ2fK/ n )( E M - l)[P(n,i;n',t-r)], (23) 

n' 

where P(n,t) is the probability of finding the system in state n at time t, and P(n, t;n',t') is the probability 
for the system to be in state n at t and in state n' at time t' . Em and Ep are raising operators acting on 
functions of n m , n p via Ejv/<?(rim, n p ) = g(n m + l,n p ) and ~Epg(n m , n p ) — g(n m ,n p + l). In the case of two-time 
quantities, e.g. P(n, t;n',t') the raising applies with respect to the second argument n'. Note that all terms 
on the right-hand side of Eq. (|23p are of order SI. This overall factor could in principle be absorbed into a 
re-scaling of time, even though we will not do so here. 

The further analysis proceeds along the lines of what was discussed for the toy model, and we will not report 
all intermediate steps in all detail. First one writes 

^rr = M(t) + tw- < 24 > 
=# = ^> + ff§, (*» 

and then systematically expands the above master equation in powers of fi -1 / 2 . To lowest order one recovers 
the mean-field equations (|16ll7p , and in next-to- leading order one finds Langevin equations for the fluctuations 
about the mean field trajectory. In the fixed-point regime of the mean-field dynamics (i.e. at large times t) 
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FIG. 5: (Colour on-line) Power spectra of the fluctuations of mRNA and protein concentrations respectively. Solid 
lines are the analytical expressions of Eqs. (|31I32[1 . Markers are from simulations at system size SI = 5000, and dashed 
lines from simulations at f2 = 500. Averages over more than 700 independent samples are taken in the simulations. 
Measurements start at t — 3000min in order to allow for equilibration. Model parameters are Po = 10, ft = 4.1, r = 
18.7, a m = ot p = 1, Hm = fi p = 0.03 @]. 



these equations read 

im{t) = -» m U(t)+a m f(P*)£ p (t-T) + ( m (t) (26) 
i P (t) = a p ^ m {t)-^ p {t)+C P (t), (27) 

where (M*,P*) is the mean-field fixed point, f(P) = df(P)/dP = -hP^il + P/P Q )- (h+1 \ ( m (t) and 

are Gaussian noise terms of zero mean and in the limit of large t, t' (when the mean-field trajectory has reached 

its fixed point) they have covariances 

(Cm(*)Cm(0> = S(t-t')[fi m M* +a m f(P*)}, (28) 

(C P (W)> = S(t-t')[fi p P* +a p M*], (29) 

(Cm(*)Cm(0> = °- (3°) 
Inverting in Fourier space one then finds after some algebraic manipulations 

£ M|2 \ = + gM(P) + fi m M*) + {a m .f'{P*)f{a p M* + y p P*) 

V "' " I (~u) 2 + n m Hp - a m a p f> '(P*) cos(o;t)) 2 + ((/x m + (i p )w + a m a p f'(P*) sin(wr)) 2 ' 

_ a 2 p (a m f(P*) + y m M *) + (uj 2 + y 2 m )(a p M* + Mp P*) 

(-w 2 + Hml^p - a m a p f'(P*) cos(wr)) 2 + ((/z m + + a m cv p /'(P*) sin(wr)) 2 ' 

The resulting power spectra for a given set of parameters used e.g. in [3| are shown in Fig. [5l and as seen 
in the figure direct simulations based on a modified Gillespie algorithm, similar to the one described in the 
section of the toy model, agree well with the theoretical predictions. Remaining discrepancies are presumably 
due to finite-size and equilibration effects. We note that the peak of the spectra shown in Fig. [5] occurs at 
an angular frequency of roughly lu 0.05 in units of 1/min , corresponding to a time period of T « 125min, 
i.e. approximately two hours and therefore close to the results from experiments reported e.g. by Hirata et 
al. [17|] . This agreement is of course a consequence of the specific choice of parameters, but it demonstrates 
that the oscillatory behaviour of mRNA and protein concentrations during somite segmentation may well be 
described as an effect of coherently amplified intrinsic noise, as opposed to cycles produced by a deterministic 
model. This enlarges the range of permissible model parameters, and our analysis as well as that of [3, [1[ 
may therefore provide additional support for the applicability of this simple stochastic model. Our theoretical 
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FIG. 6: (Colour on-line) Frequency 1/T at which the power spectrum of mRNA fluctuations has its maximum. Results 
are from Eq. (|31[1 . evaluated at the fixed point of the mean-field dynamics. The latter is obtained by integrating the 
deterministic equations ()16ll7p using an Euler-forward scheme (At = 0.1). Model parameters are Po = 100, a m = a p = 
1, jU m = Hp — 0.03. The figure is to be compared with Fig. 8 of jj]. 

analysis may also be used in order to test the robustness of the model without performing costly stochastic 
simulations throughout a large range of parameter values. In Fig. [6] for example we depict the frequency at 
which the spectrum of mRNA fluctuations has its maximum in dependence on the time delay r. This data 
is to be compared with Fig. 8 of 0], where quantitatively very similar results were obtained from actual 
stochastic simulations. Care needs to be taken though in interpreting the maximum of the power spectra as 
the frequency at which the system oscillates. The peaks in the spectra can be broad and hence several modes 
contribute. Also, of course, finite systems at small sizes may show deviations from the curves obtained from 
the system-size expansion, as the latter curves, even though they represent the next-to-leading order in f2~ 1//2 , 
can only be expected to be accurate at large system sizes. Still, Fig. El provides a theoretical confirmation of 
the findings of [J], and suggests that 2-hour cycles are found at values of the delay time of about r w 5 — lOmin 
at h — 3 and at slightly larger values of r « 15 — 17min at h = 4. It should be noted however that the variation 
of the observed time period is rather small as r and h are varied in Fig. E] so that other parameter ranges are 
not ruled out by the experimentally observed 2hr period. Still, with theoretical approaches available along 
the lines discussed in the present paper or along those of Q the most efficient way of identifying parameter 
values compatible with measurements in real-world experiments might be to use analytical expressions of the 
type presented in Eqs. (|31l32p first to narrow down the range of possible parameter values. Such closed form 
expressions can be evaluated relatively quickly and this pre-selection of model parameters based on analytical 
results is hence much less costly than performing parameter scans in stochastic simulations. Once suitable 
parameters have been identified from the theory, subsequent stochastic simulations in a much smaller range of 
parameters can then be carried out to confirm whether or not the experimentally observed behaviour is indeed 
found in the stochastic system. 

IV. DISCUSSION AND CONCLUDING REMARKS 

In summary we have successfully extended recent analyses of the effects of intrinsic noise to chemical systems 
with delay. In particular we have shown that the picture of coherent oscillations, induced by the discreteness 
of the microscopic dynamics applies to systems with delay as well. Stochastic self-sustained oscillations can 
here be found in finite systems at choices of the model parameters for which the infinite system, described by 
deterministic mean-field equations, does not exhibit cycles. This observation has important implications for the 
modelling of oscillatory biological systems with delay, as reaction rates are often not known experimentally, 
but are instead tuned in theoretical approaches, in order to ensure that simple model systems reproduce 
the experimentally observed oscillatory behaviour. Our analysis shows that confining model parameters to 
permissible ranges in which the deterministic model shows oscillations may be unnecessarily restrictive, as the 
stochastic dynamics at finite system size may well exhibit oscillatory behaviour outside these ranges of the 
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model parameters. 

While this observation as such has been made previously e.g. in 0, [E] the main contribution of the present 
work is the extension of van Kampen expansion techniques to the case of stochastic delay systems, and based 
on the resulting Langevin equation the analytical calculation of the power spectra of fluctuations about mean- 
field fixed points in delay systems. To our knowledge this has not been attempted before, even though previous 
work on Kramers-Moyal expansions in delay systems can be found in [25[ . Our approach is here complementary 
to that of who have used generating function techniques to study master equations of delay systems, and 
to computer power spectra in good agreement with simulations, but who, in our understanding, have not 
carried out a systematic expansion in the inverse system size. Since the two approaches each rely on a series 
of approximations and simplifications analytical expressions derived in the two formalisms may not be fully 
equivalent. The spectra derived in our work are however in excellent agreement with simulations, confirming 
the validity of the procedure carried out here. We have here first developed the general theory in the context 
of a simple one-dimensional system. Generalisation to more complex models with a higher number of degrees 
of freedom is possible however, and as a further example we have addressed a basic model of gene regulation. 
In particular we have studied a delay-system describing the regulatory processes underlying the expression of 
he si mRNA and Hes 1 protein. Simulational work has here recently been reported by Barrio et al. 0], and 
our work complements these mostly computational studies by an analytical computation of the spectra of the 
observed oscillations in mRNA and protein expression levels. See again also [l| for related models. Based on 
our analytical results further characterisation of the behaviour of the model is possible, without the need to 
perform computationally expensive simulations in a wide range of parameters. Our theoretical approach is 
furthermore applicable more generally, and can be expected to be useful for the theoretical understanding of 
the behaviour more intricate stochastic delay systems. 
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